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Abstract 

We analyze the statistical complexity vs. entropy plane-representation of sampled 
chaotic attractors as a function of the sampling period r. It is shown that if the 
Bandt and Pompe procedure is used to assign a probability distribution function 
(PDF) to the pertinent time series, the statistical complexity measure (SCM) at- 
tains a definite maximum for a specific sampling period £m- If the usual histogram 
approach is used instead in order to assign the PDF to the time series, the SCM 
remains almost constant at any sampling period r. The significance of tu is fur- 
ther investigated by comparing it with typical times given in the literature for the 
two main reconstruction processes: the Takens' one in a delay-time embedding, and 
the exact Nyquist-Shannon reconstruction. It is shown that £m is compatible with 
those times recommended as adequate delay ones in Takens' reconstruction. The 
reported results correspond to three representative chaotic systems having correla- 
tion dimension 2 < D2 < 3. One recent experiment confirm the analysis presented 
here. 
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1 Introduction 

The study of randomness started with Poincare, but it took almost 80 years for 
digital electronics to make enough computational resources available so that 
the investigation on the transition from order to randomness could lead to the 
fascinating issue of deterministic chaos. One hallmark of chaotic systems is the 
presence of recognizable patterns in their state space. Pattern-formation plays 
a fundamental role in the definition of i) structural complexity by Crutchfield 
[1] and ii) a number of other complexity measures, such as the so-called sta- 
tistical complexity (SCM), relevant for this paper. For an excellent review on 
complexity measures see [2]. The SCM, based on the notion of disequilibrium 
in Statistical Space, was proposed by Lopez Ruiz, Mancini and Calbet [5] 
and improved upon later in several papers |4"fo"f6 r ] . The SCM version employed 
here was advanced in In a related vein, the importance of using a "causal" 
Probability Distribution Function (PDF) to analyze the stochasticity-degree 
of chaotic and pseudo stochastic systems was emphasized in [6] and [7f8|9] . 

In this field of endeavor it is of importance to study just how the sampling 
period r-influences the final results given that digital instrumentation is widely 
used in concomitant experiments. The most common sampling criterion is 
given by the sampling theorem by Nyquist and Shannon |10fllj . It stipulates 
that the exact reconstruction of a continuous signal is possible whenever its 
Fourier spectrum vanishes for frequencies above a certain B = f ma x (the 
bandwidth). An infinite number of regularly sampled values is required and 
t must satisfy the inequality r < t NS = 1/(25) where t^s is the Nyquist- 
Shannon minimum sampling interval. The exact original continuous signal is 
recovered from the time series by using an ideal low pass filter. Note that 
in actual scenarios the number N of samples is always finite and an exact 
reconstruction is not possible. 

In a different vein, Takens [12] demonstrated that chaotic attractors generated 
by dynamical systems may be reconstructed if one has access to the infinite 
time series of regularly sampled values for just one state variable. Takens' 
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procedure allows one to obtain the main parameters characterizing the chaotic 
system, like Lyapunov exponents, dimensions, etc. For an impressive web- 
site containing very useful routines for nonlinear time series analysis see [13]. 
The ensuing reconstruction procedure uses "delayed-samples" with i T as the 
delay time, to generate a rf-dimensional vector (in a rf-dimensional embedding 
space). Note that the time delay £y must be a multiple of the sampling period 
r, since we only possess data gathered at those times. Let us remark that in 
the context of Takens' theorem the vocable reconstruction adopts a different 
meaning than that assigned by the Nyquist- Shannon theorem. In this paper we 
differentiate between a Takens reconstruction and a Nyquist reconstruction. 
We will carefully look at just how the representative point of a chaotic attractor 
in an entropy- complexity plane (H x C) changes with the sampling period. It 
will be shown that when the causal PDF prescribed by the Bandt and Pompe 
procedure [2] is used, the ensuing complexity C = C^ BP ^ exhibits a maximum 
for a given sampling period r = t^- However, such maximum does not appear 
if the conventional histogram-approach is used instead in assigning a PDF to 
the time series and thus a C = C^ hist ' is used. We call the later histogram- 
distribution a non-causal PDF. In point of fact, that of causality constitutes 
a difficult topic (see for instance [IS]). In this paper a causal PDF is one that 
takes into account the temporal correlation between successive samples, which 
entails using a PDF with a symbol assigned to each trajectory's piece of length 
L = d ■ t. 

It is of interest, for both the Takens and the Nyquist-Shannon reconstruction 
processes, to analyze the relation between i) t^ and ii) other related tempo- 
ral quantities proposed in the literature such as tr or t^s- This is one of the 
topics to be here investigated. In the case of Takens' reconstruction it is well 
known that for a finite number of samples N the quality of the reconstruction 
requires the selection of an adequate value for tr. There exist in the literature 
several different recommendations for a "good" tr, like the first zero of the 
autocorrelation function, the first minimum of the mutual information func- 
tion, etc. [16]. For an interesting discussion see [TTfTS] . In the case of Nyquist's 
reconstruction, a bandwidth criterium must be chosen in order to evaluate the 
Nyquist-Shannon time r = t NS (because the spectrum of chaotic systems is 
not band-limited). 

In Section [2] we review the Takens embedding theorem and the delay em- 
bedding reconstruction process. Section [3J deals with the Nyquist-Shannon 
reconstruction process. In Section 0] the specific SCM used in the paper is 
reviewed. Results for three paradigmatic examples: (a) the Lorenz Chaotic 
Attractor; (b) the Rossler Chaotic Attractor and (c) the chaotic attractor 
named B 7 by Chlouverakis and Sprott [19J are reported in Section [5j The first 
two attractors have been extensively studied in the literature but they both 
have a correlation dimension D2 ~ 2 while the attractor B7 has D2 ~ 2.719, 
covering almost all the 3D state-space. Conclusions and remarks are included 
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in Section [6j It is interesting to remark that the analysis performed here has 
been experimentally verified [20] . 



2 Takens' embedding theorem 

Let dx./dt = f(x) be an m-dimensional continuous dynamical system. A scalar 
measurement s(t) = s(x(t)) is a projection of the state x onto an interval 
JeR. The goal of all the embedding theorems proposed in the literature is to 
obtain a d- dimensional embedding space 21 in which x(t) can be reconstructed 
by using only s{t). This reconstruction does not need to be exact, especially in 
the case of chaotic dissipative systems that usually have attractors with box 
counting dimension Dec smaller than m, the dimension of the state space. In 
this context, reconstruction merely indicates that x(t) e 21 shares with x(t) 
some characteristics. The main requirements a reconstructed space must fulfill 
are: 

(a) Uniqueness of the dynamics in the reconstructed space. 

(b ) The reconstructed attractor must have dimensions, Lyapunov exponents 
and entropies identical to those of the original attractor. 

Consequently, the embedding of a compact smooth manifold A into M. d is 
defined to be a one-to-one map J 7 , with a Jacobian DF(x.) which has 
full rank everywhere. Let us assume that an infinite length-time scalar series 
{s n ;n = 1,2, • • • , oo} is obtained by measuring one component of the m- 
dimensional vector field at evenly spaced times t n = n-r, with r the sampling 
period. A time-delayed vector field is constructed as follows: 

™(i) = s . ™(2) ..... = m 

•^n °n j J 'n °n+l j j ^ n °n+d—l ■ V-V 

Takens proved [12] that time-delay maps of dimension 2d + 1 have the generic 
property of being the embedding of a compact manifold with dimension d, if: 
(1) the measurement function s : A — > M is (£ 2 and (2) either the dynamics 
or the measurement couples all degrees of freedom. In the original version 
by Takens, d is the integer dimension of a smooth manifold, the phase space 
containing the attractor. Thus d can be much larger than the attractor dimen- 
sion. Sauer et al. [21] were able to generalize the theorem into what they call 
the Fractal Delay Embedding Prevalence Theorem. Let Dbc be now the box 
counting dimension of the (fractal) attractor. Then, for almost every smooth 
measurement function s and any sampling time r > 0, the delay map into 
M. d with d > 2Dbc is an embedding if: (1) there are no periodic orbits of 
the system with period r or 2r and (2) there only exists a finite number of 
periodic orbits with period pr, with p > 2. Thus the main result of the embed- 
ding theorems is that it is not the dimension m of the underlying state space 
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what is important for ascertaining the minimal dimension of the embedding 
space, but only the fractal dimension Dbc of the support of the invariant 
measure generated by the dynamics in the state space. In dissipative systems 
Dbc can be much smaller than m. Let us further remark that in favorable 
cases an attractor might be reconstructed in spaces of dimension d such that 
Dbc < d < 2Dbc- For example, for the determination of the correlation di- 
mension, events of measure zero can be neglected; thus any embedding with 
a dimension larger that Dbc is sufficient. 

From a mathematical point of view, and for an infinite number of data items 
(known with infinite precision as assumed in embedding theorems), the time- 
delay tx is an arbitrary multiple of r. Thus, there exists no rigorous way of 
determining t^'s optimal value. But in a real scenario with a finite number M 
of data items, the specific value adopted by the time delay is quite important. 
Moreover, it is even unclear what properties the optimal value should have 
for best estimating a continuous system's specific property. Many different 
methods have been suggested to estimate the time-delay [T3|18II22||2"5] . In the 
case of Takens' reconstruction procedure it is possible to use a time delay 
Pr > t as pointed out above. Since in this paper the sample period r will be 
movable we will consider that both times are equal (tx = t). 

The analysis of (i) linear autocorrelations and (ii) average mutual information 
for a time series are two of the criteria often used in the literature to determine 
the best ty-region. Several characteristic times have been recommended using 
these two types of analysis. Let us recapitulate. 

(1) Time- delays induced by the discrete linear autocorrelation function. Let 
S = {s(n); n = 1, ■ • • , M} be the measured component of the vector field 
(the time series). The discrete linear autocorrelation function is a vector 
Ri defined as: 

^ M-i-l 

^ = m E \< n + «)-<*>]• [s(n)- <s>), (2) 

with < s >= J2iii s i the mean value for the time series. 

Four characteristic time delays are considered here: 

(i) The first zero crossing of Ri (if it exists). Let r and i be, respec- 
tively, the sampling period and the first zero crossing of Ri. The 
characteristic time is to = i$ r. 

(ii) The first zero crossing of R\ = Ri+i — Ri gives the first minimum 
of R. We call this characteristic time ti = i\ r. 

(Hi) The first zero crossing of R'- = R' i+1 — R\ yields the first curvature- 
change of R. We call this characteristic time t<i = ii t. 
(iv) Let % = ii/ e be the smallest i making Ri to decay to less than Rq/b. 
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The corresponding characteristic time is ty e = i\/ e r. 



(2) Delay-time induced by the discrete Mutual Information function. Let S = 
{s(n); n = 1, • • • , M} be the measured component of the vector field (the 
time series). The discrete Mutual Information function is a vector Jj, 
defined as: 



J2pki(s(n), s(n + i)) In 



k.i 



p k ,i( s ( n )^ s ( n + i )) 
Pk{s{n)) ■pi{s{n + i)) 



(3) 



Equation ([3]) is determined as follows: (a) the real interval [a, b] covered by the 
time series is partitioned into N^ ox subintervals, equal sized consecutive non 
overlapping subintervals; (b) pk is the probability to find a time series' value, 
s(n), in the k-th interval, and pk,i is the joint probability for simultaneously 
encountering a time series' value, s(n), in the k-th interval while the time 
series' value found at the 2-th posterior positions, s(n + i), falls into the the 
l-th interval. This quantity can be quite easily computed, for sufficiently small 
sizes of the partition elements (sufficiently high values of Nb ox ), and, provided 
the attractor dimension is < 2, this expression has no systematic dependence 
on N box [IH|23|. 

There exist good arguments for asserting that, if the time-delayed mutual 
information exhibits a marked minimum for a certain value of the delay ij, 
then this is a good candidate for a "reasonable" time delay. The corresponding 
characteristic time is given by tj = %i r. Nevertheless when one finds that 
the minimum of the mutual information tj lies at considerably larger times 
than ti/ e (the decay of autocorrelation function), it is worth optimizing the 
time-lag inside this range [18J. 



3 Nyquist-Shannon theorem and the minimum sampling time 



The well known Nyquist-Shannon sampling theorem [TUpTj is based on the 
Fourier Transform (FT). It states that a function s(t), t G [—00, 00], containing 
no frequencies higher than B in its FT, is completely determined by giving its 
coordinates at an infinite series of points (samples) spaced r < t NS = 1/(25) 
apart. 

Consider the case of chaotic signals that are not band-limited. A value of B 
may be defined from the Discrete FT as follows. Let {S}, be the Discrete FT 
(DFT) of {s} defined by 

M 

S k = J^ Sj e -^(i-l)(fc-l)/M (4) 
3=1 
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The power contained in the first i frequency components of the DFT's is given 
by 

PW t = £ S k SI (5) 

k=l 



where * refers to complex conjugation while the index i < M. The full power 
PW is given by Eq. fl5]) with i = M. To define B we choose the value i — i a 
in such a way as to make PWi a = aPW. The maximum frequency B for 
such a— value is then defined as B = i a /r. Once this maximum frequency is 
obtained, the Nyquist- Shannon criterion prescribes that t^g = 1/(25). Our 
present calculations were made with 0.80 < a < 0.99. 



The Statistical Complexity Measure C using Bandt and Pompe's 
prescription 



The SCM C is an informational quantifier. As such, it is a functional of a 
probability distribution function (PDF). In this case we refer to an appropriate 
PDF associated with a time series. Given the PDF, P = {pf, i — 1, • • • , iV}, 
there are several manners to obtain the functional and a full discussion of the 
subject would exceed the scope of this presentation (for a comparison amongst 
different complexity measures see the excellent paper by Wackerbauer et al. 

my 

In the present work we adopt for the SCM the functional form introduced in 
Lopez Ruiz-Mancini-Calbet seminal paper [3] with the modifications advanced 
by Lamberti et al. [5]. This functional form is an intensive (in a thermody- 
namics sense) statistical complexity C[P] given by 

C[P] = Qj[P,P e ]-H[P] , (6) 



where H denotes the amount of "disorder" given by the normalized Shannon 
entropy 

H[P] = S[P]/S maxj (7) 



and Qj is the so-called "disequilibrium", defined in terms of the extensive 
Jensen-Shannon divergence (which induces a squared metric) [5]: 

Qj[P, P e ] = Q • J[P, Pe] = Qo • {S[(P + P e )/2] - S[P]/2 - S[P e }/2} . (8) 
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In these equations P e = {1/N, - • ■ , 1/N} stands for the uniform distribution, 
while S[P] = — Ylf=i Pj hi(pj) is the Shannon entropy corresponding to the 
PDF P. We use two normalization constants, namely, 

S max = S[P e }=\nN , (9) 



and 

Qo = -2 | (^^) HN + 1) - hi(2iV) + In n\ . (10) 

The normalization constant Qo is equal to the inverse of the maximum possible 
value of J[P,P e ]. This value is obtained when one component of P, say p m , 
is equal to one and all the remaining piS are equal to zero. We have then 
< H < 1 and < Qj < 1. 

The disequilibrium Qj is an intensive thermodynamical quantity that reflects 
on the systems' architecture, being different from zero only if there exist priv- 
ileged, or more likely states among the accessible ones. C[P] quantifies the 
presence of correlational structures as well |i|5] . The opposite extremes of 
perfect order and maximal randomness possess no structure to speak of and, 
as a consequence, their C[P] = 0. In between these two special instances a 
wide range of possible degrees of physical structure exist, degrees that should 
be reflected in the features of the underlying probability distribution. 

The complexity measure constructed in this way is intensive, as many ther- 
modynamic quantities [3J. We stress the fact that the above SCM is not a 
trivial function of the entropy because it depends on two different probabili- 
ties distributions, the one associated to the system under analysis, P, and the 
uniform distribution, P e . Furthermore, it was shown that for a given H value, 
there exists a range of possible SCM values [25]. Thus, it is clear that C car- 
ries important additional information (related to the correlational structure 
between the components of the physical system) that is not contained in the 
entropic functional. 

A detailed analysis of the C-behavior demonstrates the existence of bounds to 
C that we called C max and C m i n . These bounds can be systematically evaluated 
by recourse to a careful geometric analysis performed in the space of probabil- 
ities Q [23]. The corresponding values of C max and C min depend only on the 
probability space's dimension and, of course, on the functional form adopted 
by the amount of disorder H and the disequilibrium Q. H x C diagrams are 
important and yield system's information independently of the values that 
the different control parameters may adopt. The bounds also provide us with 
relevant information that depends on the system's particular characteristics, 
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as for instance, the existence of global extrema, or the peculiarities of the 
system's configuration 

As pointed out above, the PDF P itself is not a uniquely defined object and 
several approaches have been employed in the literature to "associate" P to a 
given time series. Just to mention some frequently used P— extraction (from 
the time series) procedures, one has: a) frequency counting [26], b) time series 
histograms [S], c) binary symbolic-dynamics [27], d) Fourier analysis [28], e) 
wavelet transforms |2U|30] . f) partition entropies [51] . g) discrete entropies 
[32] . h) permutation (Bandt-Pompe) entropies [14.33J, among others. There 
is ample liberty to choose among them and the specific application must be 
carefully analyzed so as to make a good choice. Rosso et al. [6] showed that 
the last mentioned methodology may be profitably used in the plane H x C so 
as to separate and differentiate amongst stochastic, chaotic, and deterministic 
systems. It was shown in [54|5o] that temporal correlations are nicely displayed 
by the Bandt and Pompe PDF. 

Summing up, different symbolic sequences may be assigned to a given time 
series. If one symbol a of the finite alphabet 21 is assigned to each x t of the 
time series, the symbolic sequence can be regarded as a non causal coarse 
grained description of the time series because the resulting PDF will not 
have any causal information. The usual histogram-technique corresponds to 
this kind of assignment. For extracting P via an histogram one divides the 
interval [a, b] (with a and b the minimum and maximum values in the time 
series) into a finite number Nf, in of non overlapping equal sized consecutive 
subintervals : [a, b] = U^=T A i and A i = Vi ^ j. 

Note that N in equations (EJ) and (ITU]) is equal to N bin . Of course, in this ap- 
proach the temporal order of the time-series plays no role at all. The quantifiers 
obtained via the ensuing PDF-histogram are called in this paper, respectively, 
jjihist) anc j (j {hist) ^ us a j so point out that for time series with a finite al- 
phabet it is relevant to consider a judiciously chosen optimal value for A^ n 
(see e.g. De Micco et al. [8]). 

Causal information may be duly incorporated into the construction-process 
that yields P if one symbol of a finite alphabet 21 is assigned to a trajectory's 
portion, i.e., we assign "words" to each trajectory-portion. The Bandt and 
Pompe methodology for extraction of the PDF corresponds to this type of as- 
signment and the resulting probability distribution P is thus a causal coarse 
grained description of the system. Note that, in the Bandt and Pompe ap- 
proach a sort of coarse graining and word construction is effected (for method- 
ological detail see below). Note that there are other ways to get a causal PDF. 
The advantage of the Bandt and Pompe approach lies in the fact that it solves 
the problem of finding generation partitions. 
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Thus one expect that for increasing patterns' length (embedding dimension) 
the Bandt and Pompe approach retains all relevant essentials of the original 
continuous (in space) dynamics. The quantifiers obtained by appeal to this 
PDF are denoted in this paper as H^ BP ^ and C^ BP \ respectively. A notable 
Bandt and Pompe result consists in yielding a clear improvement on the qual- 
ity of Information Theory-based quantifiers |6 1 T 8 9, 34, 35, 36. 37. 38. 39, 40, 41, 42,43., 



In summarizing now the approach, note that Bandt and Pompe [T4] introduced 
a simple and robust method to evaluate the probability distribution taking 
into account the time causality of the system dynamics. They suggested that 
the symbol sequence should arise naturally from the time series, without any 
model-based assumptions. Thus, they took partitions by comparing the order 
of neighboring values rather than partitioning the amplitude into different 
levels. That is, given a time series S = {x t ;t = 1,---,M}, an embedding 
dimension d > 1 (d G N), and an embedding delay T (T £ N), the ordinal 
pattern of order d generated by 



s i — y ( x 



"s-(d-l)T) x s-{d-2)Ti - - - , Xs-T-, X s I , 



(11) 



is to be considered. To each time s we assign a d- dimensional vector that results 
from the evaluation of the time series at times s — (d — l)T, • ■ ■ , s — T, s. Clearly, 
the higher the value of d, the more information about the past is incorporated 
into the ensuing vectors. By the ordinal pattern of order d related to the time 
s we mean the permutation n = (r , r 1; • • ■ , r^-i) of (0, 1, • • • , d — 1) defined 
by 

Xs-ra^T < £s-r d _ 2 T < ' • • < Xa- n T < Xs-r T ■ (12) 



In this way the vector defined by Eq. ( II ip is converted into a unique symbol tt. 
In order to get a unique result we consider that < r»_i if x s _ ri T = ^s-r,_iT- 
This is justified if the values of x t have a continuous distribution so that equal 
values are very unusual. 

For all the d\ possible orderings (permutations) 7Tj when the embedding di- 
mension is d, their associated relative frequencies can be naturally computed 
by the number of times this particular order sequence is found in the time 
series divided by the total number of sequences, 

(t{s|s < M - d+ 1; (s) has type 7Ti} 

= W^JTi • (13) 



In the last expression the symbol Jj stands for "number". Thus, an ordinal 
pattern probability distribution P = {p(7ii),i = l,---,d!} is obtained from 
the time series. 
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It is clear that this ordinal time-series' analysis entails losing some details of 
the original amplitude-information. Nevertheless, a meaningful reduction of 
the complex systems to their basic intrinsic structure is provided. Symboliz- 
ing time series, on the basis of a comparison of consecutive points allows for 
an accurate empirical reconstruction of the underlying phase-space of chaotic 
time-series affected by weak (observational and dynamical) noise [H] . Further- 
more, the ordinal-pattern probability distribution is invariant with respect to 
nonlinear monotonous transformations. Thus, nonlinear drifts or scalings ar- 
tificially introduced by a measurement device do not modify the quantifiers' 
estimations, a relevant property for the analysis of experimental data. These 
advantages make the BP approach more convenient than conventional meth- 
ods based on range partitioning. Additional advantages of the Bandt and 
Pompe method reside in its simplicity (we need few parameters: the pattern 
length/embedding dimension d and the embedding time lag T) and the ex- 
tremely fast nature of the pertinent calculation-process [3"3"flg] . We stress that 
the Bandt and Pompe's methodology is not restricted to time series repre- 
sentative of low dimensional dynamical systems but can be applied to any 
type of time series (regular, chaotic, noisy, or reality based), with a very weak 
stationary assumption (for k = d, the probability for x t < x t+k should not 
depend on t [2]). 




The probability distribution P is obtained once we fix the embedding dimen- 
sion d and the embedding delay T. The former parameter plays an important 
role for the evaluation of the appropriate probability distribution, since d de- 
termines the number of accessible states, given by d\. Moreover, it was estab- 
lished that the length M of the time series must satisfy the condition M ^> d\ 
in order to achieve a proper differentiation between stochastic and determin- 
istic dynamics |6j . With respect to the selection of the parameters, Bandt and 
Pompe suggest in their cornerstone paper [H] to work with 3 < d < 7 with 
a time lag T = 1 . This is what we do here (in the present work d — 6 and 
T = 1 are used). Of course it is also assumed that enough data are available 
for a correct attractor-reconstruction. 

The time-causal nature of the Bandt and Pompe PDF allows for its success 
in separating chaotic from stochastic systems in different regions of the plane 

H (BP) xC (BP) p The 

main properties of the statistical complexity to be here 
employed are: (i) it is able to grasp essential details of the dynamics, because 
it employs a causal PDF, (ii) it is an intensive quantity (in the thermodynam- 
ical sense) and, (in) it is capable of discerning both among different degrees 
of periodicity and of chaos [B]. Note that H^ BF "> (the normalized Shannon per- 
mutation entropy) and C^ BP ^ (the permutation statistical complexity mea- 
sure), like other Information Theory measures (i.e. relative entropy, mutual 
information, etc.) are not system-invariant as fractal dimensions or Lyapunov 
exponents are. However, these quantifiers provide important and valuable in- 
formation on the dynamical system under analysis. For example, H^ BP ^ can 
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be consider a good approximation to the Kolmogorov-Sinai entropy [H], a 
system invariant quantity. That is, even if different time series generated by a 
nonlinear dynamical system are used for its evaluation, the global behavior of 
the system will be captured and characterized, independently of the particular 
time series chosen. 



5 Results 



5. 1 Preliminaries 



Our results refer to three paradigmatic chaotic systems with a 3-dimensional 
state space, namely, 

• The Rossler chaotic attractor, given by 



x = b + x (y — c) 
y = — x — z 
z = y + a z 



(14) 



where the parameters used here are a = 0.45, b = 2, and c = 4, correspond- 
ing to a chaotic attractor. 

The Lorenz chaotic attractor given by 



(15) 



where the pertinent parameters are a = 16, b = 4, and r = 45.92, corre- 
sponding to a chaotic attractor. 



X 


= a(y - 


x) 


< y 


= rx — 


y — xz 


z 


= xy - 


bz 



The chaotic attractor B 7 of 



x = K + z (x — ay) 
y = z [a x — e y) 



1 — x 1 — y 2 



(16) 



where K = 0.5, a = 7.0 and e = 0.23. 
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The first two systems exhibit strong differences that make them interesting 
for this study. The Rossler-oscillator "spends" most of the time near the plane 
y, z with x ~ 0. The time series obtained by sampling the x coordinate looks 
like a deJta train as Fig. [TJa illustrates. The Power Spectrum (PS) displayed in 
Fig. |2ja shows that the characteristic periodicities are very similar for all the 
three coordinates, with a higher frequency-content for x, produced by its delta 
like shape. On the other hand, the Lorenz oscillator displays a very different 
spectrum for coordinate z (see Fig. EJb), as compared with the PS of x or y. 
The reason is that the attractor spirals from the center to the border of one 
of the typical butterfly's wings and then suddenly moves to the other wing, 
as Fig. [TJb shows. Both attractors have correlation dimension D2 close to 2. 
The attractor B7 has a D2 = 2.719 ±0.156 [19], covering most of the available 
state-space. In Figs. [TJc and |2jc the time evolution of the three coordinates 
and their corresponding PS are displayed. 

The evolution of each dynamical system was determined by recourse to a 
variable-step Runge-Kutta-Fehlberg approach [49J . Evaluations were made us- 
ing sampling periods r ranging from i) 0.01 to 5, in steps At = 0.01 for the 
Rossler system, ii) 0.001 to 0.3 in steps of Ar = 0.001 for the Lorenz system, 
and Hi) 0.001 to 1 in steps of Ar = 0.001 for B7 system. With these r's one 
covers all the interesting regions: (1) the over-sampled one where r is very 
small as compared with the characteristic times considered in Sec. |2] and El 
(2) the under-sampled region where r is very large as compared with the char- 
acteristic times. For every r— value, 10 realizations were generated by starting 
from different initial conditions. The time series for the state variables x, y, 
and z were stored in a matrix of 3 x 10 5 , after skipping the first 10 4 iterations 
in order that transient states die out. 

We employed Bandt and Pompe's procedure (with dimension d = 6 and time 
lag T = 1) to assign a PDF to each time series. These PDF's exhibit an entropy 
#( BP ) and a complexity C^ BP \ Functional and C< Ms *) were evaluated 

as well, with a PDF obtained from the pertinent histogram. In this case, 
the range covered by each variable was divided into 2 16 uniformly distributed 
subintervals. The mean values < C > and < H > were computed for all the 
here considered PDF-functionals and all the realizations. For simplicity, the 
symbol < • >, meaning mean value over realizations, is suppressed. 



5.2 Our findings 

The entropy-complexity plane H^ BP ^ x C^ BP ^ is quite useful for distinguish- 
ing between stochastic noise and deterministic chaotic behavior was recently 
illustrated in [6] . In such a vein we first analyze the planar evolution of repre- 
sentative points for our systems as the sampling frequency changes (see Fig. 
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|3ja-c for the three sampled variables of the Rossler system, Fig. |3jd-f for the 
Lorenz system, and Fig. [3jg-i for the B7 system). As was explained in Sec. HJ 
"limiting curves" in the plots Fig. [3] determine the allowed values of entropy 
and complexity. There exists no PDF with entropy-complexity coordinates 
outside the region limited by these curves [25] . 

As the sample period r decreases, the representative point evolves in the 
}j( Bp ) x C^ Bps> plane, from the right planar region with low complexity and 
high entropy (under-sampling system with small sampling frequency / = 1/r) 
to the left region with low complexity and low entropy (over-sampling with 
big sampling frequency / = 1/r). The reason for such behavior is, as we 
comment in Sec. [2j that for a very big sampling period the dynamical cor- 
relations between consecutive measures are lost due to the mixing behavior 
of the chaotic system under study. That is, the measures can be completely 
uncorrelated and present characteristics like random behavior (H « 1 and 
C ~ 0). On the contrary, for very small sampling periods r the measures will 
be (strongly) time-correlated so that successive measures suggest being in the 
presence of a very ordered system (H ~ and C ~ 0). The evolution of the 
system in the H^ BP ^ x C^ BP ^ plane for regimes intermediate between these two 
extreme sampling-period instances will depend on the specific characteristics 
of the nonlinear system under study: whether it presents mixing or not; its 
chaos- degree, etc. 

It is important to keep in mind that non linear correlations are the origin of 
the geometric structures that constitute the hallmark of deterministic chaos 
[6]. Precisely, in previous works C^ BP > was seen to be a measure of the com- 
plexity induced by these non linear correlations. Our present results tell us 
that fu = - the sampling frequency corresponding to the point with 

maximum complexity - can be regarded as the optimum sampling frequency, 
i.e, the minimum sampling frequency that retains all the information con- 
cerning these structures. Higher sampling frequencies imply an oversampling, 
producing unnecessary long files to cover the full attractor, while lower sam- 
pling frequencies make the time series to lose some vital information concern- 
ing nonlinear correlations. In the oversampled case, specific <i-length ordinal 
patterns (0, 1, • • • , d — 1) and (d — 1, • • • , 1, 0) appear more frequently than 
any other ordinal pattern, causing a low entropy and a low complex- 

ity C^ BP \ Points in the region with high entropy and low complexity of the 
fj(BP) x (j(BP) p] ane correspond to undersampled trajectories of the contin- 
uous system. The samples are not correlated at all and the system behaves 
randomly. Consequently, all the ordinal patterns appear with almost the same 
frequency, as the ensuing high entropies H^ BP ^ and low complexities C^ BP ^ 
reveal. The evolution of the representative point in the plane H^ BP ^ x C^ Bps> 
near the top limiting curve is typical of chaotic systems, as demonstrated by 
the authors of 
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Let us stress the importance of using the Bandt and Pompe procedure to 
obtain the all-important PDF. A well-known result is that any chaotic map 
shares with its iterated maps the same histogram [5D]. This property is the 
basis of the skipping randomizing technique [gfoTj . For that reason De Micco 
et al. [8] proposed the use of two different PDFs, as relevant for testing (i) 
the uniformity of fi(x) (the invariant measure) and (it) the mixing constant 
"fmix pertaining to a chaotic time series, as a means to generate pseudo ran- 
dom numbers. They proceeded to employ (a) a PDF Pi based on time series' 
histograms and (b) another PDF P 2 based on ordinal patterns (permutation 
ordering) that derives from appealing to the already cited Bandt and Pompe 
method [14J. One may conjecture that continuous systems exhibit an anal- 
ogous property, viz., the sampling process does not significantly change the 
histogram of the time series and, consequently, any quantifier based on the 
PDF histogram is not sensitive to the sampling frequency. 

Such conjecture is validated in Fig. HI where the statistical complexities C^ BP ^ 
and C^ htst ^ are depicted as functions of the sampling frequency for the re- 
coordinates of the three systems (see Fig. HJa for Rossler's instance, Fig. Sib 
for Lorenz' case, and Fig. [He for B7). A similar result is obtained if C^ BP ^ and 

C (hist) 

are evaluated for the other coordinates y and z, and also if mean values 

r {BP) _ lr {BP) r {BP) r (BP) ] ,„ , r {hist) _ lr (hist) r {hist) r {hist) } ,„ 
U (xyz) — L°(a;) + °(y) + °(z) M 6 anQ U (xyz) ~ l U (x) + U (y) + U (z) I/ 6 

are used. 

Consequently, we conclude that the usual histogram technique is not useful to 
analyze the effects of the sampling frequency and that a causal PDF (here the 
Bandt and Pompe one) must be employed instead. Note that in Figs. [3] and H] 
many irregularities of the trajectories become visible. These are produced by 
the nontrivial relation between C^ BP ^ and H^ BP \ and diminish as the number 
of realizations increases. Note that in Fig. HJ where C^ BP ^ is represented as a 
function of / = 1/r, a maximum is easily detected. 

Table [T] summarizes our results, t^, t'fj , and t$ are the inverses of the sam- 
pling frequency for each coordinate producing the maximum C {BP \ We have 
evaluated t^ , tjjj, and t^ for each realization and then determined the as- 
sociated mean values and standard deviations shown in columns 1 to 6. The 
mean value of C^ BP ^ over coordinates x, y, and z also exhibits a maximum 
for the tu indicated in column 7 (< x,y,z >). All the characteristic times 
described in Sec. [2] and |3] are displayed in Table [TJ Column 8 displays the 
ratios between each characteristic time and the tu of column 7. Figures closer 
to unity in column 8 (within each class of characteristic times) are emphasized 
by using bold-types. 

As we mention in Sec. HI the normalized Shannon entropy and statistical 
complexity evaluated with PDF-Bandt and Pompe are not system invariant 
(neither are other quantifiers derived by Information Theory). The values of 
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these quantifiers depend on how the PDF is evaluated. However, they provide 
very useful information about the dynamics of the system under study. It 
is clear then that H^ BP ^ and C^ BpS) and the quantities derived from them 
could attain different values according with the specific time series considered 
(i.e. time series generated by the different coordinates of the system or a 
combination of them) even if it is reasonable to expect that their corresponding 
values do not exhibit a strong dispersion. This assumption is clearly confirmed 
if we look at the values obtained for t\j in Tabled], evaluated for different times 
series. As expected, the same kind of behavior is obtained in the case of tj 
(Average Mutual Information). 

From inspection of the results displayed in Tab. [I] we see that the time tu 
yields similar values to those obtained from the evaluation of the Average 
Mutual Information tj. The largest difference between tu and tj (see column 
8 in Tab. [T]) is observed for the B 7 system. Such disparity can be attributed 
to the fact that this system possesses an attractor with D 2 ~ 3 (that covers 
most of the state-space) in contrast with the other two systems considered 
here, the Lorenz and Rossler ones, which have D 2 ~ 2 [T8"f2"l] . 

With regards to tu and the others special characteristic times (for our three 
dynamical systems) we can state that, in the Rossler-case, such characteristic 
times are, respectively, the first minimum of the mutual information t^, the 
first zero of the autocorrelation to, and t^s, for a = 0.85. In the the Lorenz- 
instance the characteristic times closest to tu are: the first minimum of the 
mutual information i 3 , the first change of curvature of the autocorrelation t 2 , 
and t^s for en = 0.9. Finally, for B7 the characteristic times closest to tu 
are: tu e , and tpfs for a = 0.85. The first minimum of the mutual information 
appears far from the t^ —location. The shape of I m (discrete mutual informa- 
tion) is very irregular, as shown in Fig. [Hlc, which explains why such features 
takes place. 

For all the systems, the Nyquist criterion provides smaller characteristic times 
if a > 0.95 is used instead, illustrating thus on the difference between the exact 
reconstruction that is the goal of the Nyquist- Shannon Theorem and the ap- 
proximate attractor reconstruction that relies on Takens' theorem. Note that 
the autocorrelation function is based just on linear statistics and does not con- 
tain any information concerning nonlinear dynamical correlations. Kantz and 
Schreiber [18] stressed the drawback of using linear statistics for the Lorenz 
system, with a typical trajectory staying some time on one wing of the butter- 
fly attractor and spiralling then from the inside to the border before jumping 
to the other wing. The autocorrelation does not reveal the importance of the 
period of such "jumping" between wings, but the mutual information does so. 
Also, the autocorrelation of the squared signal is able to detect this periodic 
change [T8"] . 
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To assess the quality of the Takens' reconstruction process using tMi D 2 is 
evaluated for an embedding dimension ranging from d = 1 to 10 (see Fig. [7]). 
The pertinent estimation of D 2 is compared with values reported in the liter- 
ature for the continuous system [TSf52f53] . The subroutine d2 of the TISEAN 
package [13] was used for calculating D 2 . 

Let us stress that the estimation of D 2 poses a very delicate problem, so that 
we refer the readers to the documentation available for the TISEAN package 
and also to the excellent book by Kantz and Schreiber [TS]. Figs. [7J show 
that t = tM produces a time series with a D 2 that turns out to be a better 
estimation for the Correlation Dimension of the continuous system than those 
obtained with higher (under-sampled) or lower (over-sampled) values for r. 
Similar results were obtained when we consider the evaluation of the maximum 
Lyapunov exponent. 

It is interesting to remark that the time tM, for which one detects a maximum 
in the entropy-complexity plane H^ BP ^ x C^ BP \ is unique for the three dy- 
namical systems analyzed in the present work. We associate this sample-time 
with the one needed to capture all the relevant information related to the 
dynamics' nonlinear correlations. 

Recently Soriano et al. [20] theoretically and experimentally studied a chaotic 
semiconductor laser with optical feedback. This is an extremely interesting 
system because dynamical systems with time delay, like the Mackey-Glass 
or a laser with delay feedback, exhibit different relevant time scales [54"||2"0"] . 
Consequently, additional complexity-maxima could be expected (in addition 
to the one encountered in dealing with the chaotic systems of this paper). 
In the case of the laser studied in [20] the main time scales are: i) r|, the 
feedback-time providing the largest time scale, ii) Tro, that is the relaxation 
oscillation providing a shorter time scale, and finally iii) chaotic oscillations, 
that govern the fastest time scale. Soriano et al. experimentally confirmed 
that our prescription for t\j corresponds to a maximum of C^ BP \ as found 
in all the chaotic systems studied in this paper (see Fig. 10 in [20]). To our 
knowledge, this is the first controlled experiment where our prescription has 
been to hold. 

Sharing the point of view of Abarbanel [23], our prescription for the sample 
period r = t M is based on the consideration of a fundamental aspect of chaos, 
namely, the generation of information. Thus, our choice is made by taking 
into account an important property of the system we wish to describe. Re- 
mark that stable linear systems generate zero information. In consequence, 
information-generation [24] is a property of nonlinear dynamics not shared by 
linear evolutions. 

Finally, another issue to be aware of is the influence of noise. To verify the 
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robustness of the criterion here advanced, gaussian noise was added to each 
system and the value of tu was obtained as a function of the signal to noise 
S/N ratio (a) of this gaussian noise. In all cases, the value of tu remained 
constant (up to three significant decimal digits) for S/N > 20 dB. If the noise- 
strength is higher, tu suddenly decreases as the representative point in the 
plane x moves to the rightwards region (stochastic region), where 

ordering-patterns are strongly affected by noise. 



6 Conclusions 



Our main goal here was to show that a particular version of the Statistical 
Complexity Measure C^ BP \ evaluated by recourse to the probability distribu- 
tion obtained via the Bandt and Pompe procedure, allows one to determine, 
under the light of a Takens' reconstruction procedure, convenient sampling 
periods. This is done so as to get time series useful for investigating chaotic 
behavior. 

This sampling-period was called t^, and the procedure was illustrated via 
a detailed consideration of three paradigmatic systems. On the basis of these 
significant examples we conjecture that our optimality criterion may be of gen- 
eral application for chaotic systems, since C^ BP ^ is a measure of the geometric 
structures produced by nonlinear correlations, always present in this class of 
dynamical systems. We showed that tu is compatible with those specific times 
recommended in the literature as adequate delay-ones in Takens' reconstruc- 
tion. Our results closely approach the exact Nyquist-Shannon reconstruction. 
This is so because high frequencies of the Fourier Spectrum are due to chaotic 
oscillations, and this is the case for the systems studied here. 

To our knowledge, controlled experiments for different sampling times of the 
measured variables and large numbers of data have not yet been reported 
in the literature and the highest sampling frequency allowed by the digital 
acquisition system is the one usually adopted (as an exception see |20j). In 
the light of our present results it would be better, in practice, to assess the 
correct tu— value for a clever choosing of the sampling period. This permits 
one to cover the whole attractor-basin and retain in the reconstruction process 
the most important properties of chaotic systems. 

In view of the extensive use of digital acquisition systems in all kinds of exper- 
iments, we hope that experimentalists may consider the present contribution 
as a practical tool in their activities and may thus be able to validate our 
proposal. 
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Fig. 1. Evolution of the three coordinates as functions of time: (a) Rossler variables 
x, y, and z for r = 0.01 with a = 0.45, 6 = 2, and c = 4. (b) Lorenz variables x, y, 
and z for r = 0.001 with a = 16, r = 45.92, and 6 = 4. fc^ B7 variables x, j/, and z 
for r = 0.001 with K = 0.5, a = 7.0, and e = 0.23. 
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Fig. 2. Power spectra of the coordinates x, y, and z as a function of frequency: (a) 
Rossler System with r = 0.01 with a = 0.45, 6 = 2, and c = 4. for r = 0.01, and 
At = 0.01. (b) Lorenz System with a = 16, r = 45.92, and 6 = 4 for r = 0.001, and 
At = 0.001. (c) B 7 System with K = 0.5, a = 7.0, and e = 0.23 for r = 0.001, 
and At = 0.01. In all the cases the time series length considered had M = 10 5 data. 
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Fig. 3. Representation on the entropy-complexity plane: (a-c) Rossler variables x, 
y, and z for 0.01 < r < 5, At = 0.01, with a = 0.45, b = 2, c = 4. (d-f) Lorenz 
variables x, y, and z for 0.001 < r < 0.3, At = 0.001, with a = 16, r = 45.92, and 
6 = 4. (g-i) B 7 variables x, y, and z for 0.001 < r < 3, At = 0.001, with K = 0.5, 
a = 7.0, and e = 0.23. Undersampling corresponds to the high and low 

region, and oversampling to low H^ BP ' and low C^ BP \ In all the cases the 
time series length considered had M = data. The continuous lines represent, 
respectively, the maximum and minimum complexity values (for d = 6) for a fixed 
value of the entropy. 
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Fig. 3. Continuation 
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Fig. 4. Intensive statistical complexity measures C( sp ) and C( hist ^ for the coor- 
dinate x of the three systems, (a) Rossler system, (b) Lorenz system and (c) B7 
system one. C^ BP ^ is evaluated by using the PDF obtained with the Bandt and 
Pompe prescription with d = 6. C^ hlst ^ is evaluated by using the PDF obtained 
from the histogram of x, using 2 16 bins. C^ hlst>l is almost constant for all values of 
r. In all the cases the time series length considered had M = 10 5 data. The use of 
the BP prescription is of the essence in gating the complexity-maxima. The reason 
is that a change in the sampling frequency does not change the invariant measure 
of each variable of the chaotic attractor. The situation is identical to that of the 
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Fig. 5. Discrete autocorrelation function R m for: (a) Rossler System with a = 0.45, 
b = 2, c = 4. f&j Lorenz System with <r = 16, r = 45.92, and 6 = 4. (c) B7 System 
with = 0.5, a = 7.0, and e = 0.23. In all the cases the time series length 
considered had M = 10 5 data. We show both the characteristic time induced by 
C^ BP ^ (d = 6), that is tM, and its most approximate characteristic time induced by 
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Fig. 6. Discrete mutual information I rn for: (a) Rossler System with a = 0.45, 
b = 2, c = 4. (b) Lorenz System with a = 16, r = 45.92, and 6 = 4. (c) B7 System 
with K = 0.5, a = 7.0 and e = 0.23. In all the cases the time series length 
considered had M = 10 5 data. The characteristic time induced by C^ BP ^ (d = 6), 
that is tM and its closest-lying characteristic time induced by I m are shown. 



30 





Fig. 7. Mean value over realizations of the correlation dimension (< D2 >) for 
sampled x (y and z samplings produce similar results) as a function of the embedding 
dimension d: (a) Rossler System with a = 0.45, b = 2, c = 4. (b) Lorenz System 
with a = 16, r = 45.92, and 6 = 4. (c) B 7 System with K = 0.5, a = 7.0, and 
e = 0.23. In all the cases the time series length considered had M = 10 5 data. 
For oversampled series the value of D2 remains below the reference value for any 
embedding dimension. For undersampled series all linear and nonlinear correlations 
are lost and the system behaves like a stochastic one. The correlation dimension D2 
increases with the embedding dimension d. 
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Table 1 

Characteristic Times for coordinates x, y, z and their standard deviation (denoted 
by &i) of the Rossler, Lorenz and B7 chaotic attractors. System's parameters are 
a = 0.45, b = 2, c = 4 for Rossler, a = 16, r = 45.92, and b = 4 for Lorenz and K = 
0.5, a = 7.0, and e = 0.23 for B7. In all the cases the time series length considered 
had M = 10 5 data. < x,y,z > represent the corresponding time obtained averaging 
(ji BP ) over tfi e three coordinates. Each characteristic time is classified according 
with the method used to obtain it. Numbers in bold face indicate the characteristic 
time of each method that best approximates tu (the time corresponding to the 
maximum of the Statistical Complexity Measure C^ BF "> ){d = 6). 



Rossler System 


time 


X 




y 


Gy 


z 




<x,y,z> 


t/tu 


Time induced by Complexity 


tM 


1.3970 


0.0067 


1.8090 


0.1473 


1.4050 


0.0053 


1.4000 


1.000 


Times induced by Autocorrelation 


to 


1.4560 


0.0178 


1.5600 


< .00001 


1.7150 


0.0118 


1.6000 


1.1429 


ti 


3.2070 


0.0157 


3.0820 


0.0063 


3.2350 


0.0108 


3.1600 


2.2571 


*2 


0.6380 


0.0079 


1.3520 


0.0103 


1.5190 


0.0032 


1.0000 


0.7143 


h/e 


0.8480 


0.0157 


1.1700 


0.0063 


1.2930 


0.0108 


1.1200 


0.8000 


Time induced by Mutual Information 


ti 


1.6240 


0.0534 


1.4530 


0.0874 


1.5600 


0.0337 


1.5200 


1.0857 


Times induced by Nyquist-Shannon 


(0.85) 
l NS 


0.8795 




1.4645 




2.0956 




1.4799 


1.0571 


(0.90) 
l NS 


0.7292 




1.1418 




1.6661 




1.1791 


0.8422 


(0.95) 
l NS 


0.5691 




0.7447 




0.9142 




0.7427 


0.5305 
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Table 1 



Continuation. 



Lorenz System 


time 


X 




y 


Gy 


z 




<x,y,z> 


t/tM 


Time induced by Complexity 


tM 


0.0932 


0.0004 


0.0722 


0.0010 


0.0853 


0.0007 


0.0860 


1.0000 


Times induced by Autocorrelation 


to 


2.3805 


1.0482 


1.9424 


1.1457 


0.1217 


0.0018 


0.2130 


2.4767 


h 


0.4226 


0.0308 


0.3979 


0.0342 


0.2334 


0.0007 


0.2580 


3.0000 


t 2 


0.1012 


0.0006 


0.0760 


0.0001 


0.0874 


0.0005 


0.0860 


1.000 


h/e 


0.2078 


0.0308 


0.1478 


0.0342 


0.0883 


0.0007 


0.1230 


14.3023 






Time induced by Averag 


e Mutual Information 




tl 


0.1048 


0.0006 


0.0970 


0.0012 


0.0941 


0.0028 


0.0990 


1.1512 


Times induced by Nyquist-Shannon 


(0.85) 
l NS 


0.1383 




0.1082 




0.0884 




0.1116 


1.2981 


(0.90) 
l NS 


0.1138 




0.0880 




0.0670 




0.0896 


1.0419 


(0.95) 
l NS 


0.0756 




0.0639 




0.0084 




0.0493 


0.5733 
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Table 1 



Continuation. 



B7 System 


time 


X 




y 


Uy 


z 


o z 


< x,y,z> 


t/tM 


Time induced by Complexity 


tM 


0.1030 


< 0.001 


0.0930 


< 0.001 


0.9540 


< 0.001 


0.3833 


1.0000 


Times induced by Autocorrelation 


to 


0.5459 


0.3775 


0.5635 


0.4363 


1.9985 


0.0258 


1.9860 


5.1808 


h 


0.4509 


0.1332 


0.3915 


0.0957 


3.9141 


0.1161 


3.6870 


9.6182 


t 2 


0.1301 


0.0129 


0.1301 


0.0127 


1.5631 


0.0227 


0.1280 


0.3339 


h/e 


0.1813 


0.1332 


0.1747 


0.0957 


1.4504 


0.1161 


0.2830 


0.7382 


Time induced by Average Mutual Information 


ti 


1.5289 


0.2240 


1.5951 


0.2212 


1.7905 


0.1041 


1.6551 


4.3176 


Times induced by Nyquist-Shannon 


(0.85) 
l NS 


0.1901 




0.1968 




0.5319 




0.3062 


0.7987 


(0.90) 
l NS 


0.1683 




0.1736 




0.0731 




0.1383 


0.3607 


(0.95) 
l NS 


0.0539 




0.1441 




0.0076 




0.06853 


0.1787 
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